Mon. Not. R. Astron. Soc. 000, [THID (2012) Printed 6 November 2012 (MN M^jX style file v2.2) 



Carbon Cations and Silicon Atoms in the ISM: modeling 
their charge exchange reaction 
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ABSTRACT 

The time-dependent rate coefficients for the charge exchange reaction C + + Si — >• 
C + Si + for doublet and quartet states have been determined with ab initio quantum 
calculations coupled with a non-adiabatic transition model based on a simple Landau- 
Zener picture. This reaction plays a key role in determining the abundances of C, Si, 
and their ions, in the ISM since these abundances affect the fine structure cooling and 
hence the star formation rates. We also provide additional calculations to evaluate 
the differences between the gas evolution as obtained by using the empirical rate 
estimates found in the current literature and the calculations presented in this work 
which are based on our more realistic evaluation of such rates from ab initio transition 
probabilities. We shall thus show here that the new rates yield important differences 
for metal-rich environments where T < 10 4 K and the UV flux is very small, while 
becoming less important at higher T values and higher photon fluxes. 
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1 INTRODUCTION 

The chemical transformations occurring in the Interstel- 
lar Medium (ISM) and protoplanetary atmospheres plays 
a key role in explaining the evo l ution and history of its 
objects (e.g. iGalli fc Pallal ll998L iNelson fc Lariger] Il999l. 



Merlin & Chiosi 2007. 


Tesileanu et al. 20081. 


Gnedin et al. 


20091. iGlover et alJl201C 


.lYamasawa et al.ll2011 


.iGrassi et al. 



2011), since they allow to determine the total amount of 



their cooling and hence to derive specific information on 
the star formation of a given ISM region. Chemistry is also 
important for providing explanations and helping the com- 
parison betwee n the observations t hrough spectroscopi c ex- 
periments (e.g. lFuruva et ai1l2012l ; ICombes et al J 12012 ) and 
their computed simulations. Our chemical data are also em- 
ployed to compute the optical depth of the various regions 
of the ISM, thereby determining the amount of photons 
that influence the evolution of a give n environment (e.g. 
lLatif et alj[201ll ; IPetkova fc Maioll2012l) . 

Unfortunately, including an accurate and complete net- 
work of the chemical processes when endeavouring to per- 
form an ISM environment evolutionary simulation has two 
main drawbacks: the first one has a computational origin 
and it comes from the mathematical nature of the problem, 
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which is currently represented by a system of first-order dif- 
ferential equations (ODE). Solving coupled ODEs usually 
has a non- negligible computational cost for a series of rea- 
sons which are, however, beyond the aims of the present 
work so that we refer any interested reader to our recent 
study of this problem(see IGrassi et "ai] |2012l ). The second 
problem will be specifically discussed in this paper, of which 
it constitutes the main task, and is related to the accuracy 
of the rate coefficients k(T) of the molecular processes and 
chemical reactions usually employed in the model calcula- 
tions of object evolutions in the ISM. For a given reaction 
A + B — > C + D, in fact, the rate per unit volume F (i.e. 
number density per unit of time) is determined by the quan- 
tity F = k(T) riA wb, where n is the number density of the 
species indicated by the subscripts. Each reaction has there- 
fore its own k(T) that determines the efficiency of the pro- 
cess, and all the rates together affect the evolutionary model 
specifically employed. It is thus rather straightforward to un- 
derstand that to properly model the ISM one needs as many 
accurate rates as possible and that the accuracy depends on 
the method used to estimate the desired k(T). There are four 
widely-used methods to obtain a reaction rate: (i) one of the 
most accurate is to have access to experimental results, but 
unfortunately not all the reactions can be obtained trough 
this method, often because these experiments are either not 
feasible or not available at a given ti me, (ii) Another way 
is to use the Langevin model (e.g. see ISteinfeld et allll989l . 
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Sect. 8. 3) that represents a crude approximation to any re- 
action which would involve the interaction of an ion with a 
neutral species, for which is uses an adiabatic capture model 
on a single potential energy surface. It is not very accurate 
and, moreover, it misses the temperature dependence that is 
often a crucial property of ion-molecule reactions, (iii) The 
third approach is based on employing realistic chemical as- 
sumptions to describe reactions for which we have no direct 
information but are similar to others which we know more 
directly. An example would be an isotopic substitution or 
similar (for example considering that HD+H behaves the 
same as Eh+H). This method can be a necessary approxi- 
mation but it could obviously lead to large errors, (iv) The 
fourth option is the one we follow in this paper and con- 
sists in employing an ab initio method which allows us to 
obtain accurate reaction rates via quality-controlled theoret- 
ical calculations that could in principle be treated as exact 
procedures. This approach gives rate coefficients which are 
usually expected to carry fairly small errors, but it is usu- 
ally more likely to be possible for small interacting systems. 
The one we shall be discussing below constitutes a specific 
example of such a procedure. 

The process we are considering in this paper plays a key 
role when we study ISM environments with a rather weak 
local UV photon flux (e.g. where the gas is optically thick) 
and therefore the metals of such regions are not fully ion- 
ized because of the small number of ionizing photons. On 
the other hand, at the same time they are not completely 
neutral since there are other ionizing processes such as colli- 
sions and reactions with cosmic rays. For these reasons when 
the gas is optically thick we shall have to take into account 
the fccE(r) n c + nsi > where fccE(T) is the reaction rate 
coefficient for the process 



C+ + Si ^ C + Si+ . 



(1) 



The importance of determining the abundances of the dif- 
ferent metals (and ions) is related to the fact that at a 
temperature lower than 10 4 K the coolin g is dominated by 
the fine-structure cooling of the metals ijMaio et al.l 120071 : 
iGrassi et al.ll201lh . The ions of a given metal have a differ- 
ent cooli ng efficiency from that of its corresponding neutral 
partner ijSantoro fc Shulll 120061 ) . so that to employ in any 
model the correct value of £;ce(T) will affect the total cool- 
ing of the gas which is produced by that very model. In 
Sect[2] we shall discuss the computational methods which 
we employ to determine the reaction rates of the charge ex- 
change processes between C + and Si and its counterpart 
couple Si + and C. We will also provide there an analytical 
fit for their relevant kcE(T). In Sect[3]we further analyse the 
behaviour of the non-adiabatic crossings between the poten- 
tial energy curves for both quartet and doublet states of the 
combined system, we additionally present our findings for 
the cross-sections and rate constants of the charge exchange 
process over a range of temperatures and we finally com- 
pute their effects within an evolutionary model and discuss 
the differences which we find between the evolution of the 
ISM including the onl y existing empirical estimate for the 
relevant reaction rate |Le Teuff et al-lfeoOOT ) and that evolu- 
tionary behaviour obtained by employing the more realistic 
and accurate rates which have been produced in the present 
work. 



2 THEORETICAL METHOD 

The diatomic molecular cation SiC + has been the subjec t 
of ab initio studies since the work of iBruna et al] (|l98lll , 
where the dissociation channel leading to C( 3 P s )+Si + ( 2 P tl ) 
was investigated through the potential energy curves rela- 
tive to the electronic quartet and doublet states within the 
nmltireference correla tion method given by t he MRD-CI ap- 
proach. More recentlv lPramanik et al.l l|2008T ) have employed 
multireference singles and doubles configuration interaction 
to study the electronic structures and spectroscopic prop- 
erties of the SiC + cation. They have calculated all the po- 
tential energy curves leading to C( 3 Pg)+Si + ( 2 P u ) and two 
out of nine of the doublet potential energy curves which dis- 
sociate into the excited electronic state C( 1 D 9 )+Si + ( 2 P U ). 
We have used their data relative to the quartet states as a 
part of the total potential energy curves which we shall em- 
ploy below to calculate the charge transfer cross-sections 
and rate constants. The doublet potential energy curves 
for the Si + /C states have been computed within the Com- 
plete Active Space SCF ( CASSCF) methodology set up by 
IWerner fc Knowlesl l | 19851 1 The active space was given by 
nine molecular orbitals and seven electrons were kept active 
during the energy calculations. The state-averaged proce- 
dure have been used to overcame the root-flipping prob- 
lems: in particular the n-th root is calculated with a uni- 
tary weigth, and 0.01 the weigths of all the lower energy 
states. The cc-p VTZ corre l ation- consistent polarized basis 
set developed by iDunningl fl989) was employed for these 
excited state calculations and the C2v symmetry was im- 
posed to the molecular orbitals. Due to convergence prob- 
lems, the CASSCF procedure, failed in producing five out 
of the nine doublet potential energy curves associated with 
C( 1 D 9 )+Si + ( 2 P U ), and the quartet and doublet ground 
states associated with Si( 3 P 9 )+C + ( 2 P u ). A different ap- 
proach was therefore adopted for the calculations of the 
quartet and doublet potential energy curves whose asymp- 
tote is Si( 3 Pg)+C + ( 2 Pu). These last two curves have been 
calculated within the framework of the constrained Density 
Funct ional Theory (DFT) formalism of IWu &: Van Voorhisl 
(2005), using the B3LYP density functional (|Beckd Il993f ) 
and building the starting wavefunction with the positive 
charge confined to be on the carbon atom. The potential 
energy curves calculated with these two different methods 
have been asymptotically shifted in order to reproduce the 
difference in the ionization potential of the C and Si atoms, 
while the differences in the interaction regions were obtained 
from our calculations. The CASS CF calculations were per- 
formed using the MOLPRO (W erner et "all 120121 ') package, 



while the NWCHEM program (|Valiev et al 



2010) was used 



to compute the DFT potential energy curves. The range of 
the carbon-silicon distances spanned was chosen to be be- 
tween 1.25 and 6.25 A with individual steps of 0.1 A. 

The proba bility for the o ccurrence of the charge ex- 
change process l|jonhsonlll9"82T ) after every ith crossing be- 
tween the potential energy curves of different atomic charge 
localization has been calculated with in the Landau - Zener- 
Stuec kelberg (LZS) app roximation (|Landaul Il932l : IZenerl 



atueckelberg [LZit>) app 
1 19321 : IStueckelberell 19321 ) : 



of = exp 



' of 



(2) 
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where rl 



5R % 



and TQf 



7T \ Vn 



The term 



' UJ L" 1 0/l -R* 

7? :Ci is the interatomic distance at which the ith crossing oc- 
curs, v l (R Xi ) is the relative velocity between the two atoms 

at the zth crossing point. We define 8 R l — 2 " 



which depends on the derivatives of the two crossing curves 
% Vtf and 'Voo at the ith crossing point. Moreover, *Vb/ is the 
ith non-resonant exchange interaction and depends on the 
effective ionization energies, and on the ith crossing point 
coordinate: 



l V of = V fiR x 



-PR,, t (3) 

with V = P = 0.86 (-/Zsi + Vic) /V2. The 

effective ionization energies have been calculated as a func- 
tion of the effective nuclear charge Z e g as I = 0.5Z^ B . 
We have considered only the core electrons as effectively 
shielding the nucleus , and have been using the data of 
IClementi fc Raimondil (|l963h to obtain a final effective nu- 
clear charge of 3.2 and 3.3 for carbon and silicon respectively. 
The ith charge exchange probability is then employed in 
turn to calculate the corresponding ith cross-section: 

a\v) = ^Rl^ f {l~Pi )f ). (4) 

The final cross-sections are then summed over all the cross- 
ing processes and further integrated over their relevant ve- 
locities to obtain the rate constants 



k CE (T) = 4tt 



I' 



2-kKT 



v)v e 



4KT . 



(5) 



The latter quantities will then be used within the evolution- 
ary modeling described in the next Section. 

The above averaging over velocities also ensures that 
possible shifts on the crossing positions, when using different 
ab initio curves, are spread over the required velocity values 
and therefore made less sensitive in getting the final cross 
sections. 



3 RESULTS AND DISCUSSION 

3.1 Potential Energy Curves and non-adiabatic 
curve crossings 

The process of charge transfer from the cation C + to the 
neutral partner Si involves all the potential energy curves 
whose asymptotic energy levels are included between the 
two electronic states with opposite atomic charge localiza- 
tion: C( 3 P 9 )+Si+( 2 P„) and Si( 3 P s )+C+( 2 P u ). In the case 
of the quartet states all the energy levels asymptotically 
converge to either one of the two charge transfer boundary 
levels shown above (see FigfTJtop). The doublet states, on 
the contrary, present two electronic excited levels embedded 
within the two charge exchange boundary levels, namely the 
C( 1 D 9 )+Si+( 2 P U ) and C( 1 S)+Si + ( 2 P„) (see FigU bottom). 

In order to have a more complete picture of the poten- 
tial energy curves which control the charge transfer process 
we have calculated, besides the doublet and quartet elec- 
tronic ground states associated with Si( 3 P 9 )+C + ( 2 P u ), also 
the 2 2 A, the 2 2 £ + , the 3 2 A and the 2 <E> dissociating into 
C( 1 D 9 )+Si + ( 2 P U ) and the 6 2 fl and 4 2 E+ associated with 
C( 1 S)+Si+( 2 P„). We were not able, due to instabilities in 
the convergence process, to compute the other four doublet 
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Figure 1. Potential energy curves for doublet (bottom) and quar- 
tet (top) states of carbon-silicon charge exchange. The C/Si + and 
C+/Si curves have been calculated by CASSCF and constrained- 
DFT method respectively. See text for details on their calcula- 
tions. 



levels whose asymptote is C( 1 D 9 )+Si + ( 2 P U ). All the present 
calculated data turn out to be in good agreement with the 
earlier results of iPramanik et all l|2008D . 

Among the quartet states there are seven single cross- 
ings between the potential energy curve for which the posi- 
tive charge is localized on the carbon atom (black thick line 
of Fig[T]top), and the other five potential energy curves, 
while the X 4 E~ does not participate to the charge exchange 
process. The crossings are seven because the 4 II crosses three 
times the C + /Si curve, in the repulsive region at R x = 2.81 
au and in the diatomic ion bound region at R x — 3.27 au 
and R x = 3.88 au. The 4 A and the 4 E + are quasi-degenerate 
states, and cross the C + /Si curve in their repulsive parts, 
the other crossing occurs at the larger internuclear distance 
of 4.67 au. 

The doublet states crossings involve only the higher en- 
ergy curve (2 2 E~) of the six curves which dissociate into 
the ground electronic levels of the cation Si + ( 2 Pu) and neu- 
tral carbon atom ( 3 P 9 ). The crossing points, which occur at 
larger internuclear distances, involve the C + /Si curve in its 
attractive part, and the four curves, in their repulsive re- 
gions, dissociating into C( 1 D 9 )+Si + ( 2 P U ). The two higher 
energy crossing points are characterized by internuclear dis- 
tances larger than 5 au, and their energies are about 0.5 eV 
below the entrance channels C + /Si. 

The structural and energetic parameters (interatomic 
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Figure 2. Cross-sections for the charge exchange reaction in the 
doublet and quartet states of SiC + . See text for details. 



distances, energy locations and their derivatives) have been 
employed to calculate the charge exchange cross-sections a 1 
(see. Eqn. [4]) for all the crossings involving both quartet 
and doublet states. In Fig[2]we report the cross-sections as 
a function of the relative kinetic energy of the Si and C 
atoms. The behaviour of these a is almost constant below 
1 eV while at higher kinetic energies, from 1 eV to 400 eV, 
there is a marked decrease in the cross-sections with an ex- 
ponential decay factor which is almost the same for all the 
a. An exception is the behaviour of the cross-section asso- 
ciated with the quartet state II at R x = 2.81 au (bottom 
line of Fig(2|. In this case the low values of the a are due to 
the short range at which the charge exchange process can 
take place. A fact connected with the short range repulsion 
forces which therefore allow the interaction of the electrons 
with their nuclei to get stronger, thereby increasing localiza- 
tion and preventing electron charge to easily flow from one 
atom to the other. A different instance is provided by the 
crossing of 4 II at 3.27 au which produces a final cross-section 
that increases up to 10 eV and then decreases exponentially, 
similarly to the other a originating from crossings between 
the other states. This different behaviour can be explained 
in terms of the relevant potential energy curves proceed- 
ing very smoothly, and nearly parallel to each other, into 
the crossing region, a feature which shifts the maximum of 
the a to the higher kinetic energy with respect to the other 
charge exchange cross-sections between those potential en- 
ergy curves that reach more abruptly their crossing points. 

The presently computed final cross-sections have been 
employed to calculate the corresponding rate constants for 
the charge exchange process of interest, following Eqn. (0), 
as shown in FigO We also reported there the rate constants 
of the inverse charge exchange, in which the electron density 
flows from the carbon atom to the silicon cation. In Fig[3] 
we also report for comparison the rate constant which is 
usually em ployed in the ISM mo delling, taken from UMIST 
database^] (|Le Teuff et alj|200d ). We clearly see there that 
the newly computed rate constants for the direct process 
(C + +Si— >C+Si + ) are several orders of magnitude larger 
than those for the inverse process (C + +Si— »C+Si + ), and 
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Figure 3. The rate coe fficient calculated in th is Paper (solid line) 
compared to the one in iLe Teuff et~al] faOOdt) (thin dashed line). 
We also include the inverse process (thick dashed line). 



only at temperatures becoming higher than 5000 K the two 
rate constants begin to converge to very similar values. On 
the other hand , the e arlier estimate for this rate constant of 
ILe Teuff et ail (2000) does not depend on the temperature 
and therefore it is clearly higher than the rate constant cal- 
culated in this work at temperatures which are larger than 
about 1.1 x 10 4 K, while at the lower temperatures shown 
in the figure the present rate constants are more than an or- 
der of mag nitude smaller than those of the earlier empirical 
estimates l|Le Teuff et alj|2000h . 

We have also generated a fitting function that analyti- 
cally reproduce the results of our calculation as a function 
of T 

fc C E = [1.8679 x 10 s + 5.091022 x 10 10 y-°- 526561 ] - 1 (6) 
in the range T G [10, 10 5 ] K. 

3.2 Modeling the evolution of the ISM: a 
comparison 

Our newly obtained reaction rates have been included in 
a one-zone ISM model to test their direct influence over 
the corresponding evolution of the gas and relation with 
the influence already surmised in the same modelling but 
using the earlier estim ates of the same charge exchange rates 



|Le Teuff et alJl200Ch. We define fc i d as the rate coefficient 



indicated in ILe Teuff et all (|2000h . while 



fen 



1 http://www.udfa.net/ 



•Knew as the one 
provided by our calculations. In the model presented in this 
Section we follow the time-evolution of the abundances of 
the chemical species given a set of initial abundances, and a 
grid of temperature (T) and metal fraction (A) as described 
below. The chemical network is based on the more general 
one of EH over fc Jappsenl [|2007n . and we solve the system 
of differential equations ass ociated to this ne twork with a 
standard DLSODES solver l|Hindmarshlll983l ). We let the 
system evolve for 10 4 yr keeping the temperature and the 
total density both constant. 

We found that the new rates change the behaviour of 
the ISM gas only in the absence of a UV photon source, 
since the presence of such an intense source would massively 
increase the number of both Si and C ions produced by di- 
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log[T/K] log[T/K] 

Figure 4. Contour plot of logarithm of the normalized fraction defined in Eqn.lQ for low (left) and high (right) models. Contour values 
are indicated by the labels. See text for details. 



rect photoionization, thereby reducing the efficiency of the 
charge exchange due to the rapid disappearance of neutral 
partners in that reaction (i.e. Si). By keeping the latter as- 
sumption, we have therefore tested our model when select- 
ing a thick gas environment over a grid of T and A values, 
in particular by changing the gas temperature within the 
range T = [10, 10 s ] K, with different initial number densi- 
ties of the metals present in the gas according to the fol- 
lowing expression: rtx = AriHtotftx, where fix is the de- 
fault abundance of the metal X, nritot = n u + + nu + 2nn 2 
is the tot al abundance of hydroge n, while A = [0.1, 10]. 
Following Wakelam fe HerbstJ (|2008l ) we have also assumed 



"-Ho 



9 x 10" 



no 



2.56 x 10 , n c + = 1.2 x 10" 



where all the abundances are in cm 



For silicon atoms 



we have selected two different cases: n Si + = 1.7 x 10 
cm -3 and n si + = 1.7 x 10~ 5 cm" 3 labelled low and high 
model respectively. We have also chosen tih = 10 3 cm -3 , 
and n e - = Visions ni oni y a ^ the beginning of the simula- 
tions (i.e. the initial condition for free electrons). Finally, we 
have set the cosmic rays ionization rate to £cr = 1-3 x 10 -17 
s" 1 . 

To determine the differences produced within the model 
when choosing either the original rates (fcoid) or the new 
ones (fe new ) we let the system evolve for t cn d = 10 4 yr for 
both cases changing A and T within the preselected range 
indicated before, in order to fill the T— A grid of models. The 
results are plotted in Fig(4]for both low and high models as 
the logarithm of the normalized percent difference between 
the two C + abundances at t cn d, namely 



log 



'c+ 



"C+ 



(7) 



The contour lines in both panel are labelled according 
to 8 C + (e.g. S c + = — 1 represents a difference of the 10% 
between the value evolved with koid and the same for fc ne w). 
More in detail, the two plots represent in the T — A space 
the difference between the evolved n c + at t — t cn d for a grid 
of models that include fc Q id and fc now , and higher contour 



values represent higher differences, i.e. the regions of the 
T — A space where the new rate of carbon-silicon charge 
exchange calculated in this work has the larger effect. 

As expected the differences between the models essen- 
tially vanish when T « 10 4 K since at this temperature we 
have fc id ~ fc n ew ; for the same reason S c + increases at the 
lower temperatures where the two rates attain their largest 
differences. Furthermore, 5 C + is larger for the larger values 
of A, and also when the abundance of Si is larger, because 
the increase of the colliding frequencies between partners in- 
creases the probability of the carbon ion to react with the 
Si partner. 



4 CONCLUSIONS 

The work which we have discussed in the previous sections 
is directed at the determination from first principles quan- 
tum methods of the dynamics of charge transfer collisions 
between two atomic partners of the ISM composition, sili- 
con and carbon. In both cases we considered either the C or 
the Si atom to be the cationic species involved in the colli- 
sional process. The quantum mechanical calculations have 
revealed that, as expected, several ionic potential energy 
curves (PECs) are involved within the physically interesting 
range of energies and that several nonadiabatic curve cross- 
ing features are taking place, thereby guiding the charge 
transfer processes in both directions, albeit with different 
probabilities. The accuracy of the calculated PECs has been 
tested in comparison with earlier calculations and found to 
correctly describe the process of interest. Furthermore, we 
have modelled the nonadiabatic curve crossing dynamics oc- 
curring by collisional energy transfer by using the Landau- 
Zener-Stueckelberg formulae discussed in one of the preced- 
ing sections of this paper. These probabilities can further 
be used to generate Charge Exchange (CE) cross-sections 
which, in turn, yield the CE rate constants for the overall 
process and over a broad range of temperatures. The only 
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existing estimate of t he CE process has briefly presented in 
iLe Teuff et all (|200Gf i. It was essentially a Langevin-type of 
rate calculations where the atomic polarizability that drives 
the capture rate in the Langevin model (a single potential 
model) was replaced by a mass-weighted average of the two 
polarizabilities of the neutral partners in the entrance and 
exit channels of the CE process. The main feature of that 
simple model was to yield rates independent of the temper- 
ature of the ISM gas. On the contrary, the present calcula- 
tions have shown that the CE dynamics is dependent on the 
partners' velocities at the crossings and therefore generate 
by necessity a set of temperature dependent rates. We have 
further employed the new rates in a numerical modelling of 
the gas evolution, as also described in detail in the previous 
section of this work. Our results indicate that in situations 
of low-metallicity and low photon flux into the gas, the new 
rates suggest a markedly different behaviour with respect 
to the previous analysis. On the other hand, when the pho- 
ton flux increases then the overall photoionization channels 
dominate the evolution and make the collisional CE path 
essentially too inefficient. The present study has therefore 
provided a clear and specific example of the importance in 
obtaining realistic descriptions of the chemical processes ap- 
pearing in evolutionary models so that their final behaviour 
remains based on more realistic physical pictures. 
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